!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Open boundary conditions for the compressible DG code.
!!
!! \n
!!
!! The open boundary conditions are implemented here as sponge layer,
!! i.e. the flow equations are modified introducing the term
!! \f{equation}{
!!  \frac{dq}{dt} = -T(q-q_b)
!! \f}
!! where \f$q\f$ is the vector of the conservative variables,
!! \f$q_b\f$ is a boundary Dirichlet condition and \f$T\f$ is a
!! suitable damping (matrix) coefficient. Concerning the temporal
!! discretization, an implicit treatment is required, since \f$T\f$
!! can be very large. The time discretization thus reads
!! \f{equation}{
!!  (1+\Delta t\,T)q^{n+1} = q^n+\Delta t\,Tq_b
!! \f}
!! and takes the form of a weighted average between \f$q^n\f$ and
!! \f$q_b\f$. Concerning the spatial discretization, the
!! straightforward approach of discretizing the problem by
!! discontinuous finite elements turns out to be very impractical,
!! because for each element one should solve a linear system where all
!! the element degrees of freedom are coupled. In fact, the damping
!! matrix \f$T\f$ couples the variables density, energy and momentum
!! at each point, while the transformation from modal to nodal
!! representation couples the element degrees of freedom of each
!! conserved variable. To overcome this difficulty we use the
!! following algorithm:
!! <ul>
!!  <li> transform the <em>modal</em> representation
!!  \f$q^n_{\mathcal{M}}\f$ into the corresponding <em>nodal</em> one
!!  \f$q^n_{\mathcal{N}}\f$
!!  <li> solve the damping ODE in nodal space computing
!!  \f{equation}{
!!   (1+\Delta t\,T_i)\,q^{n+1}_{\mathcal{N},i} =
!!   q^n_{\mathcal{N},i}+\Delta t\,T_i\,q_{b,i},
!!  \f}
!!  where the subscript \f$i\f$ denotes a generic nodal degree of
!!  freedom (this corresponds to lumping the mass matrix in nodal
!!  space)
!!  <li> transform \f$q^{n+1}_{\mathcal{N}}\f$ to modal space
!!  \f$q^{n+1}_{\mathcal{M}}\f$.
!! </ul>
!! This algorithm reduces the computational complexity, since each
!! transformation couples either the element degrees of freedom of a
!! single conserved variable, or the pointwise values of the three
!! conserved variables. The approximation implied by this procedure,
!! equivalent to the use of a lumped mass matrix in nodal space, does
!! not represent a serious problem since accuracy is not an issue in
!! the treatment of the sponge layer.
!!
!! \note This module is controlled by the \c public variable \c
!! use_sponges: if \c coeff_sponge is associated, then \c use_sponges
!! is set to <tt>.true.</tt> and the module variables are built,
!! otherwise nothing is done. The subroutine \c sponge_apply_splitted
!! should be called <em>only</em> if <tt>use_sponges.eqv..true.</tt>.
!!
!! \todo This module relies on transformation in nodal degrees of
!! freedom, which are constructed as Lagrangian nodes. However, this
!! does not work for piecewise constant functions, since in this case
!! there are no Lagrangian nodes. One should add a specific case for
!! the piecewise constant case.
!<----------------------------------------------------------------------
module mod_dgcomp_sponges

!-----------------------------------------------------------------------

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_sympoly, only: &
   mod_sympoly_initialized, &
   ev

 use mod_linal, only: &
   mod_linal_initialized, &
   invmat

 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me, &
   t_lagnodes, lag_nodes

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_fe_spaces, only: &
   mod_fe_spaces_initialized, &
   el_nodes

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid, affmap

 use mod_dgcomp_testcases, only: &
   mod_dgcomp_testcases_initialized, &
   t_phc, phc,              &
   coeff_dir, coeff_sponge

 use mod_atm_refstate, only: &
   mod_atm_refstate_initialized, &
   t_atm_refstate,   atm_ref,   &
   t_atm_refstate_e, atm_ref_e, &
   t_atm_refstate_s, atm_ref_s

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_dgcomp_sponges_constructor, &
   mod_dgcomp_sponges_destructor,  &
   mod_dgcomp_sponges_initialized, &
   use_sponges, &
   sponge_apply_splitted

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 ! private members

! Module variables

 ! public members
 logical, protected :: use_sponges
 logical, protected ::               &
   mod_dgcomp_sponges_initialized = .false.
 ! private members
 real(wp), allocatable :: x_nod(:,:), mod2nod(:,:), nod2mod(:,:), &
   u_r(:,:,:), bbb(:,:,:,:), ccc(:,:,:,:)

 character(len=*), parameter :: &
   this_mod_name = 'mod_dgcomp_sponges'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 !> Constructor.
 !!
 !! The relevant linearized transformations for this constructor are:
 !! \f{displaymath}{
 !!  \left[ \begin{array}{c}
 !!   \Delta p \\\ \Delta T \\\ \Delta\underline{u}
 !!  \end{array}\right] = 
 !!  \left[ \begin{array}{ccc}
 !!   \frac{p_s}{\kappa}\pi_b^{\frac{1}{\gamma}} & 0 & 0 \\\
 !!   \theta_b & \pi_b & 0 \\\
 !!   0 & 0 & \mathcal{I}
 !!  \end{array}\right] 
 !!  \left[ \begin{array}{c}
 !!   \Delta \pi \\\ \Delta \theta \\\ \Delta\underline{u}
 !!  \end{array}\right]
 !!\f}
 !! and 
 !! \f{displaymath}{
 !!  \left[ \begin{array}{c}
 !!   \Delta \rho \\\ \Delta E \\\ \Delta\underline{U}
 !!  \end{array}\right] = 
 !!  \left[ \begin{array}{ccc}
 !!   \frac{1}{RT_b} & -\frac{p_b}{RT_b^2} & 0 \\\
 !!   \frac{1}{RT_b}\left(c_vT_b+\frac{u^2_b}{2}+gz\right) &
 !!   -\frac{p_b}{RT_b^2}\left( \frac{u_b^2}{2}+gz \right) &
 !!   \frac{p_b}{RT_b}\underline{u}_b\cdot \\\
 !!   \frac{1}{RT_b}\underline{u}_b & 
 !!   -\frac{p_b}{RT_b^2}\underline{u}_b & \frac{p_b}{RT_b}\mathcal{I}
 !!  \end{array}\right] 
 !!  \left[ \begin{array}{c}
 !!   \Delta p \\\ \Delta T \\\ \Delta\underline{u}
 !!  \end{array}\right],
 !!\f}
 !! from which the transformation from
 !! \f$[\Delta\pi,\Delta\theta,\Delta\underline{u}]\f$ and
 !! \f$[\Delta\rho,\Delta E,\Delta\underline{U}]\f$ is easily
 !! obtained.
 !<
 subroutine mod_dgcomp_sponges_constructor(grid,base,dt)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base
  real(wp), intent(in) :: dt

  real(wp), parameter :: tmax = 1.0e10_wp ! no need to consider larger taus
  logical :: tmask
  integer :: i, j, ie, iv, l
  type(t_lagnodes), allocatable :: lag_nod(:)
  real(wp) :: rho, uu(grid%d), u2, phi, t, pi, theta
  real(wp), dimension(2+grid%d,2+grid%d) :: ii, mtau, mma1, mma2, &
    mmap, mpam, mmm
  real(wp) :: x_nod_e(base%me%d,base%pk), tau(base%pk), taup(base%pk), &
    uuu_bnd(2+grid%d,base%pk)
  real(wp) :: gravity, rgps, frg, cv, fcv, fgamma, kappa, psk
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
     (mod_sympoly_initialized.eqv..false.) .or. &
       (mod_linal_initialized.eqv..false.) .or. &
   (mod_master_el_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) .or. &
   (mod_fe_spaces_initialized.eqv..false.) .or. &
(mod_dgcomp_testcases_initialized.eqv..false.) .or. &
(mod_atm_refstate_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_dgcomp_sponges_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   use_sponges = associated(coeff_sponge)

   if(use_sponges) then

     ! Define the Lagrangian nodes
     call lag_nodes(lag_nod,base%me%d,base%k)
     call el_nodes(base%me,lag_nod,x_nod) ! allocate x_nod
     ! notice: the number of nodes x_nod must be base%pk
     deallocate(lag_nod)

     ! Define the two transformation matrices
     allocate( mod2nod(base%pk,base%pk) , nod2mod(base%pk,base%pk) )
     do i=1,base%pk
       do j=1,base%pk
         mod2nod(i,j) = ev( base%p_s(j) , x_nod(:,i) )
       enddo
     enddo
     call invmat( mod2nod , nod2mod )

     ! Precompute the reference state at the Lagrangian nodes
     allocate( u_r(2,base%pk,grid%ne) )
     do ie=1,grid%ne
       u_r(1,:,ie) = matmul(mod2nod,atm_ref(:,ie)%rho)
       u_r(2,:,ie) = matmul(mod2nod,atm_ref(:,ie)%e)
     enddo

     ! Define the sponge coefficients
     allocate( bbb(2+grid%d,2+grid%d,base%pk,grid%ne) , &
               ccc(2+grid%d,2+grid%d,base%pk,grid%ne) )
     ! identity matrix
     ii = 0.0_wp; do iv=1,2+grid%d; ii(iv,iv) = 1.0_wp; enddo
     mtau = 0.0_wp
     ! Somehow this omp block does not work. My impression is that it
     ! is correct (in fact, there is no communication at all), but I
     ! haven't tracked this any further.
     ! !$omp parallel &
     ! !$omp   private( ie , l , x_nod_e , tmask , tau , taup , uuu_bnd , &
     ! !$omp           iv , mtau , mma1 , mma2 , mmap , mpam , mmm ,      &
     ! !$omp           rho , uu , u2 , phi , t , pi , theta , gravity ,   &
     ! !$omp           rgps , frg , cv , fcv , fgamma , kappa , psk ) &
     ! !$omp   shared( grid , base , x_nod , u_r , ii , dt , bbb , ccc ,  &
     ! !$omp           phc , coeff_sponge , coeff_dir ) &
     ! !$omp   default(none)
     gravity = phc%gravity
     rgps = phc%rgas / phc%p_s
     frg = 1.0_wp / phc%rgas
     cv = phc%cv
     fcv = 1.0_wp / cv
     fgamma = 1.0_wp / phc%gamma
     kappa = phc%kappa
     psk = phc%p_s / kappa
     ! !$omp do schedule(static)
     do ie=1,grid%ne

       ! node coordinates
       x_nod_e = affmap(grid%e(ie),x_nod)

       ! Sponge coefficients for the element
       call coeff_sponge(tmask,tau,taup,x_nod_e)

       ! Boundary state for the linearization
       uuu_bnd = coeff_dir(x_nod_e,u_r=u_r(:,:,ie))

       do l=1,base%pk ! loop on the nodes

         ! element damping coefficients: pi,theta variables
         mtau(1,1) = min( taup(l) , tmax )
         do iv=2,2+grid%d; mtau(iv,iv) = min( tau(l) , tmax ); enddo

         ! thermodynamical variables
         rho = uuu_bnd(1,l) + u_r(1,l,ie)
         do iv=1,grid%d; uu(iv) = uuu_bnd(2+iv,l)/rho; enddo
         u2 = sum(uu**2)
         phi = gravity * x_nod_e(grid%d,l)
         t = fcv * ( (uuu_bnd(2,l)+u_r(2,l,ie))/rho - u2/2.0_wp - phi)
         pi = ( rgps * rho*t )**kappa
         theta = t/pi

         ! linearized map: pi,theta -> p,T
         mma1 = 0.0_wp
         ! first line
         mma1(1,1) = psk * pi**fgamma
         ! second line
         mma1(2,1) = theta
         mma1(2,2) = pi
         ! third line (this are in fact grid%d lines)
         do iv=1,grid%d; mma1(2+iv,2+iv) = 1.0_wp; enddo

         ! linearized map: p,T -> rho,E,U
         mma2 = 0.0_wp
         ! first line
         mma2(1,1) = frg/t
         mma2(1,2) = -rho/t
         ! second line
         mma2(2,1) = frg/t*( cv*t + u2/2.0_wp + phi )
         mma2(2,2) = -rho/t*( u2/2.0_wp + phi )
         do iv=1,grid%d; mma2(2,2+iv) = rho*uu(iv); enddo
         ! third line
         do iv=1,grid%d; mma2(2+iv,1) = frg/t*uu(iv); enddo
         do iv=1,grid%d; mma2(2+iv,2) = -rho/t*uu(iv); enddo
         do iv=1,grid%d; mma2(2+iv,2+iv) = rho; enddo

         ! complete linearized map: pi,theta -> rho,E,U
         mmap = matmul( mma2 , mma1 )
         ! complete linearized map: rho,E,U -> pi,theta
         call invmat( mmap , mpam )

         ! damping matrix
         mmm = ii + dt * matmul( mmap , matmul( mtau , mpam ) )

         ! compute bbb and ccc
         call invmat( mmm , bbb(:,:,l,ie) )
         ccc(:,:,l,ie) = ii - bbb(:,:,l,ie)
       enddo
     enddo
     ! !$omp end do
     ! !$omp end parallel
   endif

   mod_dgcomp_sponges_initialized = .true.
 end subroutine mod_dgcomp_sponges_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_dgcomp_sponges_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_dgcomp_sponges_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   if(use_sponges) then
     deallocate( x_nod , mod2nod , nod2mod , u_r , bbb , ccc )
   endif

   mod_dgcomp_sponges_initialized = .false.
 end subroutine mod_dgcomp_sponges_destructor

!-----------------------------------------------------------------------
 
 !> Include the effect of the sponge layer.
 !!
 !! The sponge layer is treated implicitly in time, with an operator
 !! splitting approach with respect to the remaining tendencies.
 !!
 !! \warning In this subroutine, \c uuu does not contain the tracers.
 !!
 !! \todo One should introduce a mask to avoid looping on the elements
 !! where the damping coefficient is zero.
 !<
! pure &
 !$ subroutine omp_dummy_1(); end subroutine omp_dummy_1
 subroutine sponge_apply_splitted(uuu,grid,base)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base
  real(wp), intent(inout) :: uuu(:,:,:)
 
  integer :: ie, iv, l
  real(wp) :: uuu_nod(2+grid%d,base%pk), uuu_bnd(2+grid%d+ntrcs,base%pk)
  character(len=*), parameter :: &
    this_sub_name = 'sponge_apply_splitted'

   !$omp parallel do schedule(static) &
   !$omp   private( ie , iv , l , uuu_nod , uuu_bnd ) &
   !$omp   shared( grid , base , mod2nod , nod2mod , uuu , u_r , &
   !$omp           x_nod , bbb , ccc , coeff_dir ) &
   !$omp   default(none)
   do ie=1,grid%ne

     !------------------------------------------------------------------
     !1) Transform from modal to nodal space
     do iv=1,2+grid%d
       uuu_nod(iv,:) = matmul(mod2nod,uuu(iv,:,ie))
     enddo
     !------------------------------------------------------------------

     !------------------------------------------------------------------
     !2) Sponge layer tendency
     uuu_bnd = coeff_dir(affmap(grid%e(ie),x_nod),u_r=u_r(:,:,ie))
     do l=1,base%pk
       uuu_nod(:,l) = matmul(bbb(:,:,l,ie),uuu_nod(    :    ,l)) &
                    + matmul(ccc(:,:,l,ie),uuu_bnd(:2+grid%d,l))
     enddo
     !------------------------------------------------------------------

     !------------------------------------------------------------------
     !3) Transform from nodal to modal space
     do iv=1,2+grid%d
       uuu(iv,:,ie) = matmul(nod2mod,uuu_nod(iv,:))
     enddo
     !------------------------------------------------------------------

   enddo
   !$omp end parallel do
 
 end subroutine sponge_apply_splitted
 
!-----------------------------------------------------------------------

end module mod_dgcomp_sponges

